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Chapter  1 
INTRODUCTION 


lU _ General  Problem 

Knowledge  of  the  etrese  wave  field  created  within  the  earth 
by  an  explosive  detonation  is  required  for  the  design  of 
hardened  protective  structures  for  weapon  systems  and  for 
the  verification  of  compliance  with  nuclear  testing 
treaties.  Models  for  this  stress  wave  propagation  are 
developed  in  the  "explosive  effects  community"  by  using 
large  scale  finite  difference  codes  and  geologic  material 
models  based  on  field  and  laboratory  testing.  These  models 
are  validated  and  improved  by  comparing  calculated  stress 
and  motion  fields  with  observations  from  field  explosive 
experiments.  "Explosive  effects  community"  is  used  in  this 
study  to  represent  the  scientists  and  engineers  Involved  in 
planning,  designing,  and  analyzing  nuclear  and  conventional 
explosive  experiments. 

A  series  of  small  scale  field  experiments  conducted  over 
the  past  several  years  has  established  that  scattering  by 
random  geologic  heterogeneity  has  a  significant  influence 
on  explosively  induced  stress  wave  propagation  (Relnke  and 
Stump,  1988;  Reinke  and  Stump,  1991).  Existing  physical 
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models  of  stress  wave  propagation  in  the  explosive  effects 
community  are  based  on  the  assumption  that  the  subsurface 
geologic  material  properties  can  be  divided  into 
homogeneous  layers.  This  type  of  representation  is  often 
inadequate  as  demonstrated  by  small  soale  explosive 
experiments,  and  may  be  responsible  for  much  of  the 
discrepancy  which  remains  between  observations  from  field 
experiments  and  computational  results. 

Site  Investigation  efforts  performed  at  sites  where  field 
experiments  were  conducted  have  shown  subsurface 
heterogeneities  (Reinke  and  Stump,  1988) .  These 
heterogeneites  are  not  considered  in  the  geologic  material 
models  used  in  ground  shook  codes.  The  mean  soil 
conditions  at  test  sites  are  defined  by  constant  soil 
parameters  in  the  geologic  material  model.  While  the 
constant  soil  parameters  may  be  a  reasonable  description  of 
the  moan  conditions  at  a  tost  site,  this  simplistic 
representation  of  the  subsurface  geology  cannot  allow 
ground  shock  calculations  to  correctly  model  the 
phenomenology  which  occurs  as  the  shock  propagates  through 
the  soil  medium. 
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1-2  goal 


This  study  dsvslops  ons-dimsnsional  tschniques  for  nodsling 
ground  shock  propagation  through  spatially  random  geologic 
madia.  A  geologic  material  property  at  a  given  point  is 
unknown  until  accurately  measured,  but  it  is  not  practical 
to  perform  a  site  characterization  effort  in  enough  detail 
to  delineate  each  and  every  heterogeneity  in  the  subsurface 
soil  (Vanmarcke,  1983).  Therefore,  stochastic  modeling 
techniques  borrowed  from  the  geophysical  community  ware 
applied  to  the  problem  of  ground  shock  propagation  through 
spatially  random  geologic  media  (Frankel  and  Clayton,  1986; 
Sudicky,  1986).  Using  these  techniques,  the  spatial 
variability  (as  distinguished  from  point  variability  due  to 
imprecision  in  measurement  of  the  properties  at  a  given 
location)  in  the  material  model  is  defined  by  (1)  the  type 
of  statistical  distribution  which  defines  the  subsurface 
heterogeneity;  (2)  the  scale  or  correlation  distance  of  the 
variability;  and  (3)  the  mean  and  standard  deviation  of  the 
material  property  under  consideration.  For  the  site  used 
in  this  study,  these  statistical  parameters  have  been 
estimated  from  cone  penetrometer  testing,  laboratory 
material  property  testing,  and  seismic  surveys. 

The  site  of  Interest  for  this  study  is  McCormick  Ranch, 
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composed  of  dry  alluvium.  McCormick  Ranch  Is  a  test  site 
located  on  Klrtland  Air  Force  Base  (AFB) ,  south  of 
Albuquerque,  New  Mexico.  The  subsurface  spatial 
variability  at  McCormick  Ranch  can  be  seen  In  Figure  1. 
Figure  l  contains  a  set  6f  data  plots  from  3  cone 
penetrometer  tests  performed  on  the  sane  azimuth  (270 
degrees)  and  at  varying  distances  from  the  reference  point 
(center  of  test  bed).  Test  ov27012  was  3.66  n  (12  ft)  from 
the  tent  bed  center,  and  tests  cv27024  and  cv27040  were 
7.32  m  (24  ft)  and  12.2  m  (40  ft)  from  the  center, 
respectively.  As  shown  In  Figure  1,  a  large  degree  of 
variability  la  present  In  the  data  as  a  function  of  depth. 
The  tip  resistance  In  cv27012  varies  from  55  to  320  kg/om^ 
(782-4551  psl) ,  and  from  70  to  425  kg/om’  (996-6045  psl)  In 
CV27040.  There  la  also  a  significant  amount  of  lateral 
spatial  variability  present  as  evidenced  by  differences 
from  one  hole  to  the  next.  To  Illustrate  this,  the  tip 
resistances  from  two  arbitrary  depths  were  selected  to 
describe  the  lateral  variability: 
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ffY27012 

CV27024 
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CV27Q40 

1.22 

209 

373 

235 
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176 

180 

414 
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Depth  (m) 


Flgur*  1.  Sat  of  thraa  cona  panatromatar  soundings 
psrforasd  st  Mccornlek  Ranch,  Nsw  Naxlco. 


Figure  1  demonstrates  the  geologic  variability  present  In 
the  subsurface  material  which  Is  not  considered  In  ground 
shoclc  calculations.  The  subsurface  geologic  material 
properties  In  ground  shock  calculations  are  defined  as 
homogeneous  layers.  Figure  2  shows  a  simplified  site 
profile  of  the  McComiok  Ranch,  New  Mexico  test  site 
developed  for  use  In  ground  shock  calculations  of  high 
explosive  test  events  at  this  site  (Grant,  1988).  This 
profile  was  developed  from  a  series  of  laboratory  and  field 
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test  data  which  were  oone  penetrometer  tests,  seismic 
refraction  surveys,  and  borehole  logs.  Figure  2  shows  the 
subsurface  material  divided  Into  3  homogeneous  layers.  The 
properties  associated  with  each  layer  are  p-wave  velocity 
(Cp),  shear  wave  velocity  (C,) ,  and  soll'denslty  (p).  These 
values  are  constants,  which  do  not  account  for  the  geologic 
variability  present  in  Figure  l.  The  profile  In  Figure  2 
simply  describes  the  mean  conditions  at  the  site. 

The  main  application  of  results  from  this  study  will  be 
used  to  reduce  or  quantify  the  uncertainty  In  performing 
ground  shock/motion  calculations/predictions  based  upon 
conventional  and  nuclear  detonations. 
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Figure  2.  Slnpllfieol  site  profile  developed  for  the 
McCormick  Ranch  test  site. 


S«v«ral  studlas  hava  found  that  madivun  inhomoganalties  are 
characterized  by  correlation  analyeis.  The  application  of 
correlation  analyeis  to  the  study  of  wave  propagation  in  a 
medium  with  random  Inhomogeneities  goes  back  to  Chernov 
(1960) .  Chernov  developed  techniques  to  categorize  medium 
inhomogeneities.  He  studied  the  problem  of  amplitude  and 
phase  fluctuations  of  wave  fields  when  sound  waves 
propagated  in  the  atmosphere  and  the  sea,  the  twinkling  of 
stars,  and  the  fading  of  radio  signals. 

Correlation  analysis  has  also  been  used  to  characterize 
spatial  variation  in  earthquake  induced  ground  motions. 
The  El  Centro  differential  array  data  from  the  1979 
Imperial  /alley  earthquake  (Smith  et  al.,  1982}  ware 
analyzed  using  correlation  analysis  and  spectral  ratios  of 
array  averages.  The  differential  array  consisted  of  6  DCA- 
300  digital  recording  aocelerographs  in  a  north-south  line 
at  distances  of  0,  18.3,  54.9,  128.1,  213.4,  and  304.9  m 
(0,  60,  180,  420,  700,  and  loco  ft).  The  exact  location  of 
the  array  was  not  referenced  by  Smith  et  al.  (1982). 
However,  a  map  in  the  paper  indicated  that  the  Imperial 
fault  was  located  5  km  (3.1  miles)  from  the  first  recording 
station  (at  0  meters) .  For  the  Imperial  Valley  earthquake 


the  recording  station  at  304.9  n  (1000  ft)  did  not  trigger 
(reasons  were  not  given) .  In  the  correlation  analysis,  a 
set  of  normalized  covariance  functions  were  calculated  for 
all  possible  pairs  of  recording  stations  to  examine  the 
nature  of  the  motion  recorded  on  the  differential 
array.  The  covariance  function  measures  the  "fraction  of 
the  power  In  the  seismic  signal  that  can  be  attributed  to 
a  single  nondlsperslve  propagating  wave".  (Smith  et  al., 
1982,  p.  248).  The  degree  of  correlation  between  pairs  of 
stations  and  velocity  at  which  the  propagating  wave  moves 
across  the  array  were  obtained  from  the  covariance 
functions  of  the  accelerations.  The  normalized  covariance 
functions  for  the  P-wave  group  demonstrated  that  about  95 
percent  of  the  power  is  in  the  form  of  a  single  propagating 
wave  across  the  array  to  about  50  a  (164  ft),  but  the  power 
drops  rapidly  to  about  75  percent  at  85  m  (279  ft) .  Across 
the  entire  triggered  array,  a  distance  of  213  m  (700  ft) , 
the  normalized  covariance  Is  about  60  percent.  The  amount 
of  power  or  energy  uncorrelated  between  stations  appeared 
to  be  related  to  scattering.  Scattering  Is  caused  by  the 
Interaction  of  seismic  waves  with  the  Inhomogeneltles 
present  In  the  material  properties  of  the  geologic  medium. 

The  same  type  of  measurable  variability  was  observed  at  a 
small-aperture  seismic  array  at  Pinyon  Flat,  California 
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where  eight  local  earthquakes  were  recorded  (Vernon  et  al . , 
1991) .  The  array  consisted  of  nine  recording  stations  set 
up  in  three  nested  equilateral  triangles.  Three  stations 
were  located  300  m  (984  ft)  from  the  center  point  at  angles 
of  0,  120,  and  240  degrees.  The  middle  triangle  of 
stations  were  located  at  a  distance  of  100  m  (328  ft)  from 
the  center  point  at  60,  180,  and  300  degrees.  The  inner 
triangle  stations  were  located  30  m  (98.4  ft)  from  the 
center  point  at  30,  ISO,  and  270  degree  angles.  Pinyon 
Flat  was  selected  because  of  its  homogeneous  geology 
(granite)  and  its  planar  topography.  If  the  geology  were 
truly  homogeneous,  the  waves  would  propagate  effectively  as 
a  plane  across  the  array.  However,  this  did  not  occur. 
There  was  measurable  variability  in  the  power  spectra  and 
the  waveform  coherence  was  low  when  the  station  spacing  was 
greater  than  300  m  (984  ft)  for  wavelengths  shorter  than 
300-400  m  (984-1312  ft).  Coherence  is  a  measure  of  the 
correlation  between  two  processes  as  a  function  of 
wavelength. 

Characterizing  spatial  variability  is  also  of  great 
importance  to  groundwater  pollution  problems.  The 
complexity  of  groundwater  systems,  physical  and  chemical 
properties  of  the  subsurface  materials  that  affect  solute 
transport,  is  being  analyzed  as  spatially  random  fields. 
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Sudicky  (1986)  examined  the  spatial  variability  of 
hydraulic  conductivity  at  a  site  where  a  long  term  tracer 
test  was  performed  in  the  Borden  aquifer.  Permeability 
measurements  were  obtained  from  a  series  of  core  samples 
taken  along  two  cross  sections  at  the  Borden  tracer  test 
site.  The  samples  revealed  that  the  aquifer  was  comprised 
of  numerous  thin,  discontinuous  lenses  of  contrasting 
hydraulic  conductivity.  An  exponential  autocorrelation 
function  was  used  to  obtain  a  statistical  model  which  would 
represent  the  aquifer.  The  exponential  model  with  a 
vertical  correlation  distance  of  0.12  m  (0.4  ft),  and 
horizontal  correlation  distance  of  2.8  m  (9.2  ft)  closely 
approximated  the  estimated  autocorrelation  functions.  The 
statistical  description  of  the  Borden  aquifer  was  used  to 
determine  the  dispersion  of  the  Injected  tracer,  and  then 
compared  to  the  dispersion  rates  observed  during  the  long¬ 
term  tracer  test.  Because  the  results  were  consistent, 
Sudicky  believes  that  a  statistical  understanding  of  random 
hydraulic  parameters  will  provide  meaningful  estimates  of 
transport  parameters  in  other  aquifers. 

Correlation  functions  were  used  by  Frankel  and  Clayton 
(1986)  to  study  the  scattering  of  elastic  and  acoustic 
waves  in  two-dimensional  media  with  random  spatial 
variations  in  seismic  velocity.  The  paper  was  concerned 
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with  variations  In  seismic  velocity  In  the  earth's  crust 
with  scale  lengths  (correlation  distances)  ranging  from 
tens  of  kilometers  down  to  tens  of  meters.  The  scattering 
model  consisted  of  a  random  velocity  perturbation  applied 
to  a  homogeneous  velocity  field.  The  perturbation  fields 
(also  called  random  media  in  the  paper)  were  characterized 
by  a  correlation  function »  correlation  distance,  and 
standard  deviation.  The  Gaussian,  exponential,  and  Von 
Karman  correlation  functions  were  used  to  construct  the 
random  media. 

The  method  used  by  Frankel  and  Clayton  to  construct  the 
random  media  Is  of  great  Importance  to  this  study  because 
a  similar  procedure  was  utilized.  The  procedure  used  to 
construct  the  random  media  on  a  two-dimensional  grid  was 
outlined  in  the  paper  as  follows:  "First,  a  random  number 
generator  assigned  a  velocity  v(x,z)  sequentially  to  each 
point  of  the  grid.  The  random  velocity  field  was  then 
Fourier  transformed  to  wave  number  space,  filtered  to 
achieve  the  desired  spectrum,  and  transformed  back  to  the 
spatial  domain  to  yield  the  velocity  field  for  the 
simulations"  (Frankel  and  Clayton,  1986,  p.  6469) .  The 
results  from  the  simulations  were  then  compared  to  actual 
observations  to  constrain  models  of  the  crustal 
heterogeneity.  The  results  demonstrated  that  the  Gaussian, 
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exponential,  and  Von  Kaman  random  media  with  correlation 
distances  greater  than  10  km  (6.2  miles)  produce 
correlation  functions  similar  to  the  actual  observations. 

This  literature  review  concentrated  on  the  usage  of 
correlation  analyses  in  various  disciplines.  The 
correlation  analysis  was  used  to  characterize  the  spatial 
variation  in  strong  ground  motions,  examine  the  spatial 
variability  of  hydraulic  conductivity,  and  study  the 
scattering  of  elastic  and  acoustic  waves  in  media  with 
random  spatial  variations  in  seismic  velocity. 

JLti _ 0Ja.1ftatiygg.,.g£ 

To  achieve  the  goal,  the  specific  objectives  of  this  study 
were  to  do  the  following: 

1.  Develop  a  one-dimensional  random  geology 
generator  to  generate  the  random  geologic  variability 
profiles  to  be  used  in  finite  difference  calculations. 

2.  Modify  a  ground  shock  finite  difference  code  to 
incorporate  the  random  geologic  variability  factors. 
The  factors  will  be  used  to  perturb  the  average  soil 
material  properties  to  generate  random  geologic 
material  property  profiles  that  represent  the 
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Inhomogeneltles  present  in  the  subsurface  material. 

3 .  Perform  linear  elastic  calculations  using  the 
modified  code. 

4.  Compare  the  results  of  the  calculations  with 
closed  form  solutions  and  with  experimental  data. 


i±s _ SM.cigiqg  .Qg,  ThU  stutiiY 

This  research  is  unique  within  the  explosive  effects 
community  in  the  following  ways.  First,  stochastic 
techniques  are  used  to  produce  spatial  perturbations  In  the 
material  model.  These  perturbations  simulate  the 
subsurface  heterogeneities.  Existing  solutions  for  plane 
wave  propagation  In  a  heterogeneous  material  have  been 
obtained  by  using  probabilistic  techniques  (Rohani,  1982), 
which  do  not  incorporate  the  concept  of  spatial 
variability.  Second,  the  spatial  random  geologic 
variability  factors  are  directly  incorporated  into  the 
AFTON-ID  finite  difference  ground  shock  code  which  has  not 
been  done  previously  (Trulio,  1966). 
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The  organization  of  this  thesis  Is  as  follows:  chapter  1 
Is  the  Introduction  which  states  the  problem,  goal,  review 
of  literature,  and  objectives.  Chapter  2  Is  Theory  and 
Methodology.  The  topics  covered  In  chapter  2  are 
autocorrelation  functions,  random  geology  generator, 
theoretical  and  computed  autocorrelations,  AFTON-ID  finite 
difference  code,  and  code  modifications  and  calculations. 
The  results  and  analysis  are  presented  in  chapter  3.  In 
chapter  4,  conclusions  are  dra%m  based  on  the  results  from 
the  various  calculations  and  recommendations  are  given  for 
future  work  in  this  area. 
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Chapter  2 

THEORY  AND  METHODOLOGY 

2.1  Autocorralatlon  Function 

Random  or  atoohaatlo  prooaaaaa  ara  commonly  characterized 
through  autocorrelation  functions  or  spectral  density 
functions.  The  spectral  density  function  Is  simply  the 
Fourier  transform  of  the  autocorrelation  function  of  a 
random  process. 

The  autocorrelation  function  Is  defined  as  a  measure  of  the 
degree  of  linear  relationships  between  neighboring  values 
In  a  random  process  (Bendat  and  Plersol,  1971) .  The 
general  autocorrelation  fimctlon  Is  defined  as 

t+t)  -  -i  jgXjf(t)Xjf(t+T)  (a»l) 

N  -  number  of  samples 

where,  R^x  ■■  autocorrelation  function  of  x^Ct) 
t  reference  point 
r  -  spatial  separation 

The  quantity  R„  Is  always  a  real-*valued  even  function  with 
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a  maximum  at  r«0,  and  may  be  either  positive  or  negative. 
In  terms  of  the  autocorrelation  function,  the  mean  value  of 
x(t)  is  given  by 


(2.2) 


Another  form  of  the  autocorrelation  function  la  the 
normalized  autocorrelation  function  (correlation 
coefficient) 

^VSFxlJWySFlejlJlE^ 

-1  i  R^(x)  i  1  (2.4) 

When  the  autocorrelation  function  is  at  a  pealc,  the  random 
process  has  a  high  degree  of  linear  correlation  with  itself 
at  t  and  t-fr.  When  the  autocorrelation  function  is  near 
zero,  the  random  process  is  uncorrelated  with  itself  at  t 
and  t-fr. 

The  autocorrelation  function  is  used  in  many  disciplines, 
which  include  geotechnical  engineering,  statistics, 
geophysics,  and  hydrology.  The  autocorre?atlon  function  is 
commonly  referred  to  as  simply  the  correlation  function  in 
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the  literature.  Throughout  this  paper,  the  terms 
autocorrelation  and  correlation  will  be  used 
Interchangeably . 

The  most  common  autocorrelation  functions  used  to 
characterize  random  geologic  variability  In  the 
selsmologlcal  and  hydrological  communities  are  the 
Gaussian,  exponential,  and  von  Karman  functions  (Frankel 
and  Clayton,  1986;  sudlcky,  1986).  The  correlation 
function  and  corresponding  l^D  power  spectra  equations  are 
shown  In  Table  l.  The  correlation  functions  In  Table  1  are 
theoretical  eqfuatlons  for  a  Gaussian,  exponential,  and  von 
Karman  distribution,  respectively,  and  are  always  positive. 

The  data  set  from  McCormick  Ranch  used  to  statistically 
characterize  the  subsurface  heterogeneity  suggests  that  the 
exponential  and  Gaussian  correlation  functions  are  most 
appropriate  for  describing  the  spatial  variability. 
Therefore,  this  study  has  focussed  on  the  Gaussian  and 
exponential  correlation  functions  to  define  the  subsurface 
heterogeneities.  The  correlation  functions  and  1-0  power 
spectra  are  plotted  In  Figure  3.  The  1-D  power  spectra 
plots  were  generated  with  a  correlation  distance  (a)  of  50 
cm. 
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For  Gaussian  and  exponential  correlation  functions  shown  In 
Table  1,  the  correlation  distance  (a)  marks  the  lag  (r) 
whore  the  correlation  function  has  the  value  of  e‘‘. 
Therefore,  when  r"a  the  normalized  lag  (r/a)  Is  equal  to 
R^>iQ.37  as  shown  In  Figure  3. 


Table  1.  Correlation  functions  and  1-D  power  spectra 
equations. 


Correlation 

Function 

1-D  Power 

Spectra 

GAUSSIAN 

e-xV«‘ 

av/i:  «•**•*/< 

EXPONENTIAL 

0-x/a 

2a 

l+k*a* 

von  KARMAN 

K^U/a) 

a 

(l+k®a*) 

Where , 

r  ■ 

offset  (or  spatial  lag) 

a  ■ 

correlation  distance 

k  - 

wave 

number,  ko27r/\,  where 

X  Is  wavelength 

Ko  ■■  Bessel  function 
e  ">  exponential  function 
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AUTOCORRELATION  FUNCTIONS 


NORMALIZED  LAG  (r/a) 

1-D  POWER  SPECTRA 


10“^  10"^  10"^  10° 
IfAVE  NUMBER  (k) 


Flgur*  3.  Gaussian  and  Exponsntlal  thsorstloal 
autocorrslatlon  functions  and  1-D  powsr  spectra. 
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A  tttchniqua  was  developed  to  construct  spatially  random 
geologic  variability  factors  which  would  be  Incorporated 
into  the  AFTON-10  finite  difference  code  (Trulio,  1966)  to 
simulate  the  subsurface  inhomogene  it  lea.  The  technique  has 
been  named  the  random  geology  generator.  The  procedure 
used  to  develop  the  random  geology  generator  is  as  follows. 

1.  First,  the  desired  correlation  function,  f,  is 
constructed  In  the  wave  number  domain  using  the 
equations  (l-D  power  spectrum)  in  Table  1. 

2.  Next,  a  random  sequence,  y,,  is  generated  using  a 
random  number  generator  (discussed  later  in  this 
chapter) .  This  random  sequence  is  then  transformed  to 
the  wave  number  domain  by  means  of  a  fast  Fourier 
transform. 

3.  A  random  phase  term,  9,  is  obtained  by 

tf-tan'‘[ijaag(V(k))/roaI(V(k))  ],  where  Y(k)  is  the 
Fourier  transform  of  the  sequence  y|.  0  is  known  as 

the  phase  angle  of  the  Fourier  transform,  and  imag  and 
real  stand  for  the  imaginary  and  real  part  of  the 
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Fourier  transform,  respectively. 


4.  Since  a  complex  number  (x+ly)  can  be  written  as 


X  *  ly  •  \]x^  *  a  * 
the  random  sequence  Is  written  as 


(a. 5) 


a^*  <a.«) 

where  X(k)  Is  a  random  sequence  In  the  wave  number 
domain  having  the  desired  correlation  function  (f ) . 
Transforming  X(k)  back  to  the  space  domain  yields  X|, 
which  Is  a  random  sequence  In  space. 


5.  The  random  sequence  Is  then  standardized  to  a 
specif io  mean  and  standard  deviation.  The 
transforming  equation  is 

(^- (»•’) 

where  and  are  the  mean  and  standard  deviation  of 
the  random  sequence  X|,  and  are  the  desired  mean 
and  standard  deviation,  and  x.  Is  the  sequence  of 
random  geologic  variability  factors. 


The  program  to  generate  the  random  geologic  variability 
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factors  Is  listed  in  Appendix  A.  The  input  parameters 
needed  to  run  the  program  are  1}  the  desired  correlation 
function  (Gaussian  or  exponential),  2)  number  of  random 
variables  to  be  generated,  3)  the  desired  mean  and  standard 
deviation,  4)  correlation  distance,  and  5)  a  seed  for  the 
random  number  generator.  The  correlation  distance  in  the 
program  is  defined  in  terms  of  the  number  of  sample 
intervals. 


2.2.1  correlation  Distance 

The  correlation  distanoe  is  a  stochastic  parameter  used  to 
measure  the  size  of  variations  in  a  data  set.  The 
correlation  distance  is  defined  as  a  measure  of  the 
separation  distance  where  two  measurements  are  essentially 
independent  of  each  other. 

The  correlation  distanoe  used  in  this  study  is  based  on 
work  done  by  Jewell  (1988).  A  cone  penetrometer  data  set 
from  the  McCormick  Ranch  test  site  on  Kirtland  AFB  composed 
of  alluvium  was  analyzed  for  stochastic  and  statistical 
properties.  Various  autocorrelation  functions  were  used  to 
determine  the  spatial  variability  in  the  subsurface 
properties,  correlation  distances  of  30-50  cm  provided  a 
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good  fit  between  the  model  and  data.  For  this  study,  two 
correlation  distances  of  30  and  50  cm  were  Investigated. 


2.2^2 _ Random  Variability  Factors 

The  Gaussian  and  exponential  functions  wore  used  to 
generate  the  random  variability  factors.  The  variability 
factors  are  used  to  perturb  the  average  soil  material  modal 
to  generate  random  geologic  material  property  profiles  that 
represent  the  Inhomogeneltles  in  the  subsurface  material  In 
the  ground  shock  code.  Since  the  grid  in  the  ground  shock 
code  has  a  limit  of  soo  cells  or  zones  and  the  number  of 
variability  factors  generated  must  equal  2*  for  the  random 
geology  generator,  n  was  chosen  to  be  9.  Therefore,  512 
variability  factors  were  generated. 

Figure  4  shows  the  variability  factors  constructed  with  the 
Gaussian  function.  The  variability  factors  constructed 
with  the  exponential  function  are  shown  in  Figure  5.  The 
profiles  shown  in  Figures  4  and  5  represent  the  subsurface 
haterogeneltlas  with  depth.  The  variability  factor 
profiles  were  computed  with  different  zone  sizes. 
Therefore,  in  order  to  compare  the  profiles  directly,  the 
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VARIABILITY  FACTORS  VARIABILITY  FACTORS 


GAUSSIAN  FUNCTION 
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Flgur*  4.  Variability  factor  profilaa  for  tha  Gauaaian 
function  at  dlffarant  zona  slzas,  and  with  raapaot  to 
nomalizad  dlatanoa. 
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^ARIABILmr  FACTORS  VARIABILITY  FACTORS 


EXPONENTIAL  FUNCTION 
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Flgur*  5.  Variability  factor  profllaa  for  the 
Exponantlal  function  at  dlffarant  zon*  alzaa,  and  with 
raspact  to  normal Izad  dlatanca. 
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variability  factors  are  shown  with  respect  to  the 


nonnalized  distance  in  the  figures.  The  same  input 
parameters  were  used  to  generate  the  different  profiles. 
When  these  generated  profiles  are  compared  with  the  cone 
data  profiles  in  Figure  1,  they  appear  quite  realistic  (see 
Figure  6) . 

The  variability  factors  from  the  exponential  function  show 
more  high  frequency  roughness  or  noise  than  those  from  the 
Gaussian  function.  The  roughness  is  due  to  the  fall  off 
rate  in  the  spectral  amplitude  between  the  two  correlations 
(see  Figure  3).  Both  spectra  are  flat  out  to  a  specific 
wave  number  (k) ,  but  at  higher  wave  numbers  the  exponential 
falls  off  at  where  N  is  the  number  of  space 
dimensions  (Fisk  et  al.,  1991).  The  amount  of  roughness  in 
the  profile  is  controlled  by  the  fall  off  rata  of  the 
spectra.  The  spectrum  with  more  energy  at  higher  wave 
numbers  shows  more  roughness  in  the  generated  variability 
profiles. 


_ Random  Number  Generator 

The  random  number  generator  from  the  MATLAB  software 
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Figure  6.  Comparison  between  cone  penetrometer  tests  and 
generated  subsurface  variability  simulations. 
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package  (PC-MATLAB/  1989)  was  utilized  in  this  research  (in 
Section  2.2).  MATLAB  is  a  high-performance  interactive 
software  package  used  for  scientific  and  engineering 
numerical  computations.  MATLAB  provides  an  easy-to-use 
environment  where  problems  and  solutions  are  expressed 
almost  exactly  as  they  are  written  mathematically.  A 
specific  programming  language  is  not  required.  Typical 
uses  for  MATLAB  are  general  purpose  numerical  computations, 
algorithm  prototyping,  and  solving  problems  with  matrix 
formulations  that  arise  in  disciplines  such  as  automatic 
control  theory,  statistics,  and  digital  signal  processing. 

The  MATLAB  function  RAND  Is  a  uniform  random  number 
generator  based  on  linear  oongruential  methods.  The 
algorithm  for  this  random  number  generator  was  proposed  by 
D.  H.  Lehmer  (1951).  The  algorithm  Involves  the  choice  of 
two  fixed  integer  parameters 

(1)  modulus:  a  -  a  large  prime  integer 

(il)  multiplier:  a  -  an  integer  in  the  range  2, 

3,....,m  —  1 

and  the  generation  of  the  integer  sequence  Zi,  Z^,  Z^,...,Z^.^ 
is  defined  by  the  recursive  formula 

(ill)  Z|  "•  (aZi.|)  (mod  m)  for  i  <■  1,  2,... 
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where  Zq  Is  the  seed  or  starting  value.  To  obtain  the 
desired  random  numbers,  U,  on  (0,1),  the  sequence  of  Z's  Is 
normalized  by 

(Iv)  U|  -  Z|/m  for  1  1,  2,  ... 

The  values  of  a  and  a  must  be  selected  In  the  manner  that 
the  function,  f(Z),  la  a  full  period  generating  function, 
and  the  full  period  sequence,  Zi,  Z2,...,Z^,,  Is  random.  The 
modulus  a  ■  2^'  *•  1  was  suggested  by  Lehmer  (1951).  The 
multiplier  a  <■  7^  «  16807  was  suggested  by  Lewis,  Goodman 
and  Miller  (1969),  based  on  the  fact  that 

f(Z)  -  168072  mod  2147483647  (8*5) 

Is  a  full  period  generating  function.  This  function  has 
also  demonstrated  evidence  of  randomness  through  various 
testa  (Lehmer,  1951). 

Random  number  generators  can  be  subjected  to  empirical  and 
theoretical  tests  to  determine  how  well  the  generated  Uj's 
resemble  values  of  a  true  Independent  and  Identically 
distributed  (IID)  random  variable.  In  addition  to  the 
statistical  tests  conducted  In  by  Lehmer  (1951),  the  chi- 
square  (empirical)  test  of  uniformity  was  applied  to 
equation  (2.5)  In  this  study.  The  remits  from  the  chi- 
square  test  indicate  the  U/s  generated  do  not  behave  In  a 
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way  which  is  significantly  different  from  what -would  be 
expected  from  truly  independent  Identically  distributed 
random  variables.  The  computations  and  results  from  the 
chi-square  test  are  included  in  Appendix  B. 


Tha  Gaussian  and  exponential  correlation  functions  were 
used  to  generate  the  random  geologic  variability  factors. 
The  random  factors  are  applied  to  the  homogeneous  soil 
material  model  to  develop  synthetic  heterogeneous  soil 
profiles.  In  order  to  determine  the  validity  of  the  random 
geology  generator,  the  autocorrelations  of  the  variability 
factors  were  compared  to  the  theoretical  autocorrelations. 
The  theoretical  autocorrelation  functions  shown  in  Table  1 
were  developed  baaed  on  integral  Fourier  transform  theory 
from  -00  to  +«  over  space.  In  order  to  construct  a 
numerical  realization  (random  geologic  variability 
factors) ,  a  discrete  process  was  used.  The  integral  is  no 
longer  from  -«o  to  +oe,  but  rather  from  0  to  512.  Therefore, 
the  adequacy  of  the  computed  autocorrelations  obtained  from 
the  discrete,  limited  realizations  was  compared  to  the 
theoretical  autocorrelations. 
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Figure  7  shows  ths  thsoratlcal  and  computed  auto¬ 
correlations  for  the  Gaussian  and  exponential  functions, 
respectively,  for  a  correlation  distance  of  SO  cm.  Figure 
8  shows  the  autocorrelations  for  a  correlation  distance  of 
30  cm.  The  characteristics  associated  with  the  random 
factors  In  Figures  7  and  8  are  a  mean  of  l,  a  standard 
deviation  of  10  percent,  and  zone  size  of  10  cm.  The 
random  number  generator  seed  was  sat  at  100. 

The  comparison  between  the  theoretical  and  computed 
autocorrelations  In  both  eases  was  extremely  good. 
Therefore,  one  would  conclude  that  the  random  geology 
generator  constructing  the  discrete,  limited  realizations 
Is  working  correctly. 


The  AFTON  1-D  (one  dimensional)  finite  difference  code  was 
used  In  this  study.  The  AFTON  l-D  code  Is  a  finite 
difference  ground  shock  code  designed  to  solve  problems  In 
the  field  of  dynamic  contlnutim  mechanics  (Schuster  el  at., 
1984} .  AFTON  1-D  may  be  used  in  any  of  the  three  1-D 
symmetric  modes;  plane,  cylindrical  or  spherical.  Figure 
9  shows  the  three  symmetric  modes.  The  only  direction  of 
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EXPONENTIAL  FUNCTION 


NORMAUZED  LAG  (r/a) 

Figure  7.  Theoretical  and  computed  autocorrelations  for 
the  Gaueelan  and  Exponential  functions  for  a  correlation 
distance  of  50  cm. 
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GAUSSIAN  FUNCTION 


Correlation  Distance  =  30  cm 


NORMAUZED  LAG  (r/a) 


EXPONENTIAL  FUNCTION 


NORMALIZED  LAG  (r/a) 

Figure  8.  Theoretical  and  conputed  autocorrelatlone  for 
the  Gauealan  and  Exponential  functions  for  a  correlation 
distance  of  30  cm. 
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Flgur*  9.  ThrM  l-D  aynunstric  modaa  for  tha  AFTON-ID 
coda.  (From  Sohustar  at  aX.,  1984) 
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motion  allowed  in  the  thraa  modas  Is  that  psrpendlcular  to 
ths  shaded  area  (Aj)  In  the  figure.  The  media  are  divided 
Into  discrete  cells,  Identified  by  J,  and  the  code 
cycles  through  the  equation  of  motion  for  all  active  J's 
using  the  output  of  the  previous  cycle  to  advance  each 
tlmestep.  The  state  of  stress  In  each  cell  or  zone  was 
calculated  by  a  linear  elastic  model.  The  parameters  that 
define  the  linear  elastic  model  are  material  density  (p  - 
1643  kg/m^) ,  compresslonal  wave  velocity  (o  <■  508  m/s) ,  and 
Poisson's  ratio  (i*  «  0.25). 

This  code  vras  used  to  perform  homogeneous  and  random 
calculations  in  the  spherical  mode  with  a  pressure  step 
function  (p  ■  1  MPa)  applied  at  the  grid  boundary 
(spherical  cavity  wall)  to  simulate  ground  shock 
propagating  through  the  grid  from  an  explosive  source.  K 
zone  size  of  10  cm  was  used  throughout  thn  grid. 


2«iLJL _ validatrian,  of.  gfllgulfttiQnajL.  Ttghniaui 

As  a  check  on  thb  validity  of  the  c^lculatlonal  technique, 
the  homogeneous  finite  difference  calculations  were 
compared  to  a  closed  form  solution  of  wave  motion 
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davelopad  by  Sharp*  (1942) .  Sharpe  state*  the  problem  as 
follows:  "Given  a  spherical  cavity  of  radius  a  within  a 
homogeneous,  ideally  elastic,  infinite  medium  of  density  p 
and  oompressional  wav*  velocity  o;  to  find  the  elastic  wave 
motion  which  results  from  application  of  an  arbitrary 
pressure  P(t)  to  the  Interior  surface  of  the  cavity" 
(Sharpe,  1942,  p.  146).  The  problem  is  Illustrated  In 
Figure  10. 

The  primary  reason  for  the  waveform  eomparlsona  between  the 
homogeneous  finite  difference  and  closed  form  calculations 
was  to  validate  both  the  parameter  setup  of  tha  code,  as 
well  as  the  geometric  setup  of  tha  finite  difference  grid. 


Figure  10.  Illustration  of  the  idealized  :!iodel. 
(From  Sharps,  1942,  p.  146} 
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A  computer  proqram  was  written  for  Sharpe's  closed  form 
solution  and  Is  Included  In  Appendix  C  along  with  the 
equation  derivations.  The  program  was  written  to  generate 
velocity  and  displacement  time  histories  at  several 
distances  from  the  center  of  the  cavity. 

Comparisons  between  the  closed  form  solution  and  the  AFTON 
calculations  are  shown  In  Figures  11  through  13  at 
arbitrary  ranges  of  200,  500,  and  900  cm  from  the  cavity 
center,  respectively.  The  vertical  scales  In  the  figures 
vary  to  show  the  maximum  velocities  and  displacements  at 
the  different  ranges.  Also  shown  in  Figures  11  through  13 
are  the  time  histories  generated  with  different  zone  sizes. 
In  order  to  oheoK  the  influence  of  frequency  and  boundary 
effects  on  the  calculations  due  to  the  zone  size,  several 
zone  sizes  were  compared.  The  homogeneous  calculations 
were  performed  with  2,  5,  10,  and  20  cm  zone  sizes. 

The  velocity  and  displacement  waveforms  from  the  closed 
form  solution  and  AFTOK  calculations  are  similar.  The 
waveform  agreement  increased  as  the  distance  from  the 
cavity  wall  also  increased.  The  reason  for  the  agreement 
is  explained  by  the  simplified  displacement  equation  used 
by  Sharpe  (1942) .  The  simplified  equation  is  adequate  for 
displacements  at  distances  more  than  a  few  times  the  radius 
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Flgura  11.  Comparisons  bstwssn  ths  olossd  form  solution 
and  AFTON-ID  calculations  for  2,  5,  10,  and  20  cm  zonss 
at  200  om  from  ths  cavity  osntsr. 
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Flgura  12.  Conparlsons  batwaan  tha  cloaad  tom  solution 
and  AFTON-ID  calculations  for  2,  5,  10,  and  20  cm  zonts 
at  500  cm  from  tha  cavity  cantar. 
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Figure  13.  Conparlsons  between  the  closed  form  solution 
and  AFTON-ID  calculations  for  2,  5,  10,  and  20  cm  zones 
at  900  cm  from  the  cavity  center. 
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of  the  cavity.  The  cavity  radius  for  the  calculations  was 
set  at  175  cm.  Based  on  the  agreement  between  the  closed 
form  solution  and  AFTON  calculations,  one  would  conclude 
the  parameter  setup  used  In  the  code  Is  adequate. 

As  shown  In  Figures  11  through  13,  the  difference  In  zone 
size,  between  the  2,  s,  and  10  cm  zones.  In  the  finite 
difference  grid  was  not  significant.  However,  numerical 
noise  was  present  In  the  20  cm  rone  calculation,  and  the 
velocity  peak  values  were  lower.  Therefore,  a  zone  size 
greater  than  10  cm  would  not  be  recommended.  Based  on 
these  comparisons,  the  geometric  setup  of  10  cm  zones  for 
the  homogeneous  and  random  calculations  was  used  throughout 
this  study. 


2u.S _ AFTQN  l-D  Code  Modifications 

In  order  to  determine  the  Influence  of  spatial  geologic 
variability  on  ground  motion,  the  AFTON  l-D  finite 
difference  code  was  modified  to  Incorporate  the  random 
geologic  variability  factors.  The  variability  factors  were 
assigned  to  each  zone  or  cexl  on  the  grid.  They  were  then 
used  to  perturb  the  average  soil  material  properties  to 
represent  the  inhomogeneltles  In  the  subsurface  material. 
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Th«  AFTON  1-D  code  consists  of  53  subroutines  which  include 
functions  for  initializing  constants,  setting  up  the 
problem,  adding  new  cells  as  the  problem  progresses, 
setting  the  new  tlmestep,  and  storing  output  data.  The 
subroutine  of  interest  for  this  study  was  SUBROUTINE 
ESINIT.  The  parameters  from  the  linear  elastic  model  are 
read  in  at  this  subroutine  and  assigned  to  the  zones.  This 
subroutine  was  modified  so  that  the  set  of  variability 
factors  could  be  read  in  to  the  material  model  and  applied 
to  the  parameters  in  each  zone.  Since  only  a  linear 
elastic  material  model  is  being  explored  in  this  study,  the 
material  parameters  perturbed  by  the  variability  factors 
are  density,  Poisson ^s  ratio,  and  seismic  wave  speed.  The 
material  parameters  were  assigned  to  500  zones  and  then 
varied  by  SOO  variability  factors  composing  a  Gaussian  or 
exponential  distribution. 
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Chapter  3 

RESULTS  AND  ANALYSIS 


hi _ Random  Calculatlona 

The  random  calculations  vara  ganaratad  with  tha  modified 
AFTON  l-D  finite  difference  coda.  In  tha  first  phase  of 
tha  analysis,  tha  statistical  properties  of  tha  variability 
factors  ware  varied  to  examine  tha  affects  on  tha  computed 
output,  velocity  and  displacement  time  histories. 
Correlation  distances  of  30  and  50  cm  ware  compared,  and 
the  standard  deviation  In  tha  variability  factors  varied 
from  0  to  20  percent. 


IxIjlI _ Correlation  Distance  comparisons 

The  variability  factors  used  In  the  random  calculations  to 
compare  the  difference  In  correlation  distance  were 
generated  with  an  arbitrary  standard  deviation  of  10 
percent,  and  Initial  seed  of  100.  The  correlation 
distances  of  30  and  50  cm  are  compared  in  Figures  14 
through  17.  The  waveforms  in  Figures  14  through  17 
correspond  to  ranges  of  200,  500,  900,  and  2000  cm  from  the 
center  of  the  cavity,  respectively. 
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Flgura  14.  Corralatlon  dlatano*  comparisons  for  30  and 
50  cm  at  a  rangs  of  200  cm. 
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Flgur*  15.  Corralation  dlatanca  conparlaona  for  30  and 
SO  cm  at  a  rang#  of  500  cm. 
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Figure  16.  correlation  distance  comparlaone  for  30  and 
50  ont  at  a  range  of  900  cm. 
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Flgura  17.  correlation  dlstanoa  comparisons  for  30  and 
50  cm  at  a  range  of  2000  cm. 
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The  wavefonts  generated  with  the  two  correlation  distances 
vary  slightly,  but  there  are  no  significant  differences. 
The  time  of  arrivals  are  the  same  for  both  cases  at  the 
different  ranges. 


3.1.2  standard  Deviation  and  Initial  Seed  Comparisons 

The  differences  in  standard  deviation  and  initial  seed  in 
the  random  calculations  are  shown  in  Figures  IB  through  21. 
The  three  waveforms  in  each  plot  represent  the  soil 
properties  modeled  as  a  homogeneous  half  space  (sd-0) ,  and 
as  random  geologic  material  property  profiles  generated 
with  a  50  cm  correlation  distance,  different  random  seeds 
(s-tlOO,  sm456)  ,  and  standard  deviations  (sd«5%,  sd"10%, 
sd>15%,  and  ad'"20%} .  These  waveforms  were  generated  at  a 
range  of  500  cm  from  the  center  of  the  cavity.  Figures  18 
through  21  represent  5,  10,  15,  and  20  percent  standard 
deviations,  respectively. 

The  velocity  time  history  waveforms  generated  from  the 
homogeneous  calculations  (sd<«0}  show  the  initial  peak  and 
then  decay  to  zero.  For  the  time  histories  generated  with 
the  random  variability  profile  calculations,  the  waveforms 
do  not  decay  as  rapidly  with  time,  but  rather  show  the  late 
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Flgur*  18.  Random  calculations  gsnaratad  with  5  parcant 
standard  davlatlon  in  variability  factors. 
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Flgur*  19.  Random  calculatlona  ganaratad  with  10  parcant 
atandard  davlation  in  variability  factors. 
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Figur*  20.  Random  calculations  gsnsratsd  with  15  percent 
standard  deviation  in  variability  factors. 
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Flgur*  21.  Random  calculations  ganaratad  with  20  parcant 
standard  deviation  In  variability  factors. 
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time  noise  that  is  often  present  In  measured  field 
measurements.  The  other  significant  difference  between  the 
homogeneous  and  random  calculations  In  Figures  18  through 
21  is  the  reduction  in  the  peak  velocity  value.  As  the 
level  of  variability  in  the  soil  properties  Increases  (from 
0  to  20  percent  standard  deviation) ,  the  attenuation  of  the 
ground  shock  motion  increases. 

Also  Illustrated  in  Figures  18  through  21  are  the 
differences  in  the  random  calculations  duo  to  the  initial 
seeds  used  in  the  random  number  generator.  The  random 
variabi]ity  factors  generated  by  the  two  initial  seeds 
produced  significant  differences  in  the  velocity  and 
displacement  waveforms.  Even  though  the  two  sets  of  random 
variability  factors  are  statistically  the  same,  the  initial 
seed  used  in  the  random  number  generator  influences  the 
computed  velocity  and  displacement  waveforms.  In  order  to 
account  for  the  differences  produced  by  the  initial  seeds, 
representative  random  waveforms  (mean  and  mean  ±  1  standard 
deviation)  were  generated  based  on  work  performed  by 
Vanmarcke  (1979) . 

Vanmaroke  developed  a  methodology  to  probabilistically 
characterize  the  variability  of  geotechnical  parameters  and 
design  sampling  and  testing  programs  to  maximize  the 
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probabilistic  information.  Based  on  Vanmarcke's  findings 
and  experi'^nce,  he  suggests  that  a  minimum  of  6-8  data 
points  are  needed  to  estimate  means  and  variances  or 
coefficient  of  variation  with  a  reasonable  degree  of 
reliability.  Therefore/  representative  random  waveforms 
were  generated  using  lo  different  seeds. 


In  order  to  compare  the  results  produced  by  the 
Introduction  of  spatial  subsurface  variability  to 
homogeneous  ground  shock  calculations  and  experimental 
data/  representative  random  waveforms  (mean  and  mean  1 
standard  deviation)  were  generated.  The  representative 
random  waveforms  were  generated  with  a  50  cm  correlation 
distance,.,  a  Gaussian  distribution/  soil  properties  \iith  a 
10  percent  standard  deviation/  and  by  using  3  sets  of  10 
different  seeds  (Set-l/  Set-2/  and  Set-3}.  The  seeds  used 
in  the  3  sets  are  shown  In  Table  2. 

r.epresentativa  velocity  waveforms  were  generated  at  varying 
ranges  from  the  pressure  cavity  center  (200,  500,  900/  and 


Table  2.  Three  sets  of  random  number  generator  seeds. 


234 

1000 

572 

345 

2000 

802 

456 

3000 

33 

567 

4000 

534 

678 

5000 

499 

789 

6000 

955 

890 

7000 

748 

980 

8000 

554 

1234 

9000 

624 

2345 

10000 

89 

2000  cm)  for  the  3  sets  of  independent  seeds.  The 
waveforms  are  presented  in  Figures  22  through  27. 

The  soil  properties  in  these  random  calculations  were 
varied  by  10  percent  based  on  down  hole  seismic  cone 
investigations  conducted  at  the  McCormick  Ranch  test  site 
by  The  Earth  Technology  Corporation  (Final  Report  to  Air 
Force  Weapons  Laboratory,  Contract  No.  F8PR0171120007, 
1987)  . 


The  representative  random  velocity  waveforms  presented  In 
Figures  22  through  27  are  shown  In  Figures  28  and  29  as 
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time: 


TIME  (M8«e](MlO*) 


Figure  22.  R«pr«s«ntatlv«  valoclty  wavaforms  from  Set-l 
at  200  and  500  cm  from  tha  cavity  cantar. 
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Figure  23.  R«prMantatlv«  valoclty  wavaforns  from  Sat-l 
at  900  and  2000  cm  from  tha  cavity  cantar. 


58 


VELXITY  Ct»/psacl  t*10  )  VELOCITY  (co/iisec)  I  *10 


Figurs  24.  Rttprassntatlv*  vvloolty  wavaforms  from  Sat-2 
at  200  and  500  cm  from  tha  cavity  cantar. 
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Figure  25.  Rapresantatlvs  valocity  wavafoma  from  Sat-2 
at  900  and  2000  cm  from  tha  cavity  cantar. 
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Figur*  26. 

Raprasantativa  valoclty  wavaforms  from  Sat-O 

at  200  and 

500  cm  from  tha  cavity  cantar. 
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Figure  27.  Representative  velocity  waveforms  form  set-s 
at  900  and  2000  cm  from  the  cavity  canter. 
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RANGE  (cm) 


Flgura  28.  Random  and  honoganaoua  paak  valoclty 
amplitudaa  varaua  ranga  for  Sat~l. 


Flgura  29.  Random  and  hom^^d^a^ui^^^k  velocity 
anplltudea  versus  range  for  Set-'2  and  Set-3, 
respectively. 
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peak  velocity  amplitudes  versus  range  for  Set-1,  Set-2,  and 
Set-3,  respectively.  This  type  of  data  presentation  is 
useful  for  analyzing  the  decay  rate  of  the  random  and 
homogeneous  calculations  for  the  three  cases.  Figures  28 
and  29  show  the  moan,  mean  -t-l  standard  deviation,  mean  -1 
standard  deviation,  and  homogeneous  peak  velocity 
amplitudes. 

In  general,  the  peak  velocity  amplitudes  from  the 
homogeneous  calculation  are  larger  than  the  mean  amplitudes 
and  in  some  cases  are  not  within  the  ±l  standard  deviation 
bounds.  These  plots  indicate  that  the  ground  shook 
attenuation  rate  is  larger  for  the  random  calculations 
relative  to  the  homogeneous  oaloulations. 


3.1 2.1,2 _ Coefficients  of  Variation 

The  peak  velocities  obtained  from  the  random  and 
homogeneous  calculations  were  normalized  to  the  mean 
amplitudes  for  the  three  random  cases  to  eliminate  a  biased 
comparison  between  the  sets  ot  independent  seeds.  The 
ratio  is  known  as  the  coefficient  of  variation  (Bendat  and 
Piersol,  1971). 
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The  coefficients  of  variation  for  the  peak  velocity 
amplitudes  are  shown  in  Figure  30  as  a  function  of  range. 
The  mean  +i  standard  deviation,  mean  -1  standard  deviation, 
and  homogeneous  peak  velocity  amplitudes  were  normalized  to 
the  mean  amplitude  at  the  different  ranges  for  each  set. 
The  mean  ±1  standard  deviation  for  the  three  sets  are 
represented  by  the  open  triangle,  diamond,  and  circle 
symbols.  The  homogeneous  velocity  amplitudes  normalized  to 
the  mean  of  each  set  are  represented  by  the  closed 
triangle,  diamond,  and  circle  symbcls  for  set>l,  set-2,  and 
Set-3 ,  respectively . 

The  widths  of  the  ±1  standard  deviation  bounds  or  error 
bounds  are  consistent  for  the  three  sets  of  seeds,  except 
for  Set-1  at  the  200  cm  range.  The  error  bound  for  Set-1 
is  approximately  23  percent  from  the  mean,  where  Set-2  and 
and  Set-3  are  only  13  percent  from  the  mean.  However,  the 
widths  of  the  error  bounds  at  the  different  ranges  do  vary. 
The  widths  of  the  largest  error  bounds  at  the  different 
ranges  are  23,  18,  25,  and  14  percent  for  the  200,  500, 
900,  and  2000  cm  ranges,  respectively.  The  average  error 
bound  is  20  percent  for  all  ranges. 

Figure  30  also  shows  the  homogeneous  peak  velocity 
amplitudes  with  respect  to  the  mean  amplitude  and  bounds. 
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A  Set  1  ^  Set  2  o  Set  3 
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Figure  30.  Coeffloiente  of  variations  with  respect  to 
range  for  Set-1,  Set-2,  and  Set-3.  The  open  syidDols 
represent  the  bounds  for  the  3  sets  of  random 
oaleulations,  and  the  closed  symbols  represent  the 
homogeneous  calculations  norma] ized  to  the  mean  of  each 
set. 
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The  hoTnogenenous  amplitudes  for  the  three  cases  are  closest 
to  the  mean  and  within  the  bounds  at  the  200  cm  range.  At 
the  500  and  900  cm  ranges,  the  homogeneous  amplitudes 
normalized  to  the  Set>3  means  are  outside  the  bounds.  At 
the  2000  cm  range,  the  homogeneous  amplitudes  for  the  three 
oases  are  outside  the  bounds.  The  average  homogeneous 
amplitude  Is  30  percent  from  the  mean,  where  the  largest 
error  bound  is  14  percent  at  2000  cm. 

Figure  30  shows  that  the  ground  shock  decay  rate  with  range 
for  the  homogeneous  caloulatlons  is  smaller  than  the  random 
calculations. 


3.2.3  other  Random  Waveform  Effects 

In  the  previous  sections,  only  the  peak  amplitudes  from  the 
homogeneous  calculations  were  compared  to  those  of  the 
representative  random  calculations.  In  this  section,  the 
complete  homogeneous  velocity  waveforms  will  be  compared  to 
the  mean  random  velocity  waveforms. 

Four  homogeneous  velocity  waveforms  are  shown  in  Figure  31 
as  a  function  of  range.  The  mean  random  waveforms  for  Set- 
1,  Set-2,  and  Set-3  are  shown  in  Figures  32  through  34, 
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Figure  32.  Representative  mean  velocity  waveforms  as  a 
function  of  range  for  Set-1. 


70 


Figur*  33.  Rttprasantativ*  aaan  valoclty  wavaforaa  as  a 
function  of  ranga  for  Sat-2. 
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Figure  34.  Representative  mean  velocity  waveforms  as  a 
function  of  range  for  Set-3. 
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respectively.-  The  mean  random  waveforms  which  take  into 
account  the  spatial  subsurface  variability  show  a  reduction 
in  the  peak  amplitudes  as  well  as  a  shift  to  a  lower 
frequency.  The  energy  is  spread  out  with  time  in  the 
initial  positive  cycle  for  the  mean  random  waveforms.  The 
mean  random  waveforms  also  show  the  late  time  ringing  or 
noise  often  seen  in  experimental  data. 


3.3  Experimental  Data  Comparison 

In  the  final  phase  of  the  analysis,  the  random  calculations 
were  compared  to  experimental  data.  The  analysis  was  more 
qualitative  than  quantitative  for  the  following  reasons: 

1.  The  random  calculations  performed  in  this  study 
were  purely  linear  elastic  calculations.  All  near 
source  explosion  data  involve  non-linearity. 

2.  The  calculations  were  1-dimensional,  where  the 
actual  conditions  are  3 -dimensional. 

When  an  explosive  test  is  buried  sufficiently  deep  (no 
surface  effects) ,  the  experiment  can  be  assumed  to  be  l- 
dimensional.  Therefore,  the  random  calculations  were 
compared  in  a  qualitative  way  to  three  approximately  1- 
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dimensional  explosive  tests.  These  three  tests  are  MERLIN, 
MP-2,  and  CHEAT. 


3.3.1  MERLIN  Test 

MERLIN  was  a  contained  underground  nuclear  experiment 
conducted  in  desert  alluvium  in  Yucca  Flats,  Nevada  Test 
Site  (Ferret  and  Bass,  1975).  MERLIN  was  a  10  kiloton  test 
at  a  depth  of  296  m  (971  ft)  with  numerous  gages  placed  at 
shot  depth  approximating  a  one-dimensional  configuration. 

The  instriimentation  layout  for  the  MERLIN  test  is  shown  in 
Figure  35.  The  ranges  at  stations  U3,  U4,  U5,  and  U7  are 
350,  700,  1100,  and  2503  ft,  respectively,  from  the  charge. 
The  velocity  waveforms  are  shown  in  Figure  36.  The 
velocity  waveform  at  U7  was  obtained  from  an  accelerometer 
measurement  and  not  a  velocity  measurement,  therefore,  the 
symbol  J  was  placed  in  the  plot.  As  seen  in  Figure  36,  the 
late  time  ringing  or  noise  introduced  into  the  random 
calculation  by  the  subsurface  variability  (see  Figures  32- 
34)  is  also  seen  in  the  MERLIN  results  for  the  various 
ranges.  In  previous  analyses  of  data  of  this  typo,  the 
late  time  ringing  has  been  attributed  to  reflections  from 
various  interfaces,  or  from  surface  spall  impact  (Murphy, 
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KPIHIO 


HORIZONTAL  RANGE  (K) 

Figure  35.  MERLIN  teet  inetrujnentation  layout. 
(From  Parrat  and  Baaa,  1975) 


75 
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TIME  (sec; 

Flgur*  36.  MERLIN  results  at  ranges  350,  700,  1100,  2503 
ft  from  the  charge  (top  to  bottom) . 

(From  Ferret  and  Bass,  1975) 
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1991} .  The  random  calculations,  however,  demonstrate  that 
a  relatively  small  amount  of  spatial  variability  in  the 
material  properties  can  produce  similar  effects. 


3.3.2  Material  Properties  Test 

The  MP-2  test  consisted  of  a  20  ton  nitromethane  charge 
buried  20  m  below  the  surface  in  dry  alluvium  at  the  Yuma 
Test  Site,  Yuma,  Arizona  (Rinehart  and  Veyera,  1986) .  A 
total  of  178  gages  were  placed  in  this  test  on  5  different 
azimuths.  The  results  from  this  test  were  used  to  study 
the  variation  in  peak  particle  velocities  at  approximately 
the  same  ranges  for  5  azimuths. 

Figure  37  shows  two  velocity  plots  with  5  waveforms  for 
approximate  ranges  of  9.5  m  (top  plot)  and  12.6  m  (bottom 
plot) .  These  plots  show  the  variability  in  recorded 
measurements  at  the  same  distance,  but  different  azimuths. 
The  peak  particle  velocities  at  the  9.5  m  range  vary  from 
20.5  to  27.5  m/sec,  with  an  average  velocity  of  23.4  m/sec 
and  a  coefficient  of  variation  of  11  percent.  The  peak 
particle  velocities  at  the  12.6  m  range  vary  from  4  to  18 
m/sec,  with  an  average  velocity  of  11.6  m/sec,  and  a 
coefficient  of  variation  of  14  percent.  Those  ±1  standard 
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rigur«  37.  velocity  tin#  hi»tori««  from  MP-2  at 
approxinataly  9.5  m  (top)  and  12.6  m  (bottom)  ranges. 
(From  Rinehart  and  Veyera,  1986) 
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deviation  bounds  are  similar  to  the  uncertainty  bounds  for 
the  random  calculations  shown  In  Figure  30. 


3.3.3  contained  Hloh  Explosive  Alluvium  Test  f CHEAT) 

The  CHEAT  experiment  consisted  of  a  253  lb  TNT  charge 
detonated  at  a  depth  of  11.5  m  In  dry  alluvium  at  McCormick 
Ranch  (Stump  and  Re  Inks,  1987}  at  a  site  near  where  the 
cone  penetrometer  study  was  conducted.  The  Instrumentation 
layout  Is  shown  In  Figure  38  (top  plot)  for  the  CHEAT 
experiment . 

The  only  measurement  of  Interest  for  this  study  Is  located 
at  a  depth  of  7.5  m  from  the  surface,  since  all  the  other 
measurements  are  close  to  the  free  surface.  The 
acceleration,  velocity,  and  displacement  time  histories  are 
shown  In  Figure  38  (bottom  plots) .  The  velocity  time 
history  for  this  test  again  shows  the  late  time  noise  or 
ringing  seen  In  the  random  velocity  waveforms.  Stump  and 
Relnka  attributed  some  late  time  features  of  the  waveforms 
to  free  surface  effects,  but  again  the  random  calculations 
show  that  a  small  amount  of  spatial  variability  can 
Introduce  similar  late  time  features  Into  the  waveforms. 
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Flgum  38.  CHEAT  instrumentation  layout  (top)  and  tima 
hiatoriaa  (bottom) .  (From  Stump  and  Rainka,  1987) 
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Chapter  4 
CONCLUSIONS 


One-dlmenglonal  techniques  for  modeling  ground  shock, 
propagation  through  spatially  random  geologic  media  were 
developed  in  this  study.  The  subsurface  spatial 
variability  was  defined  by  a  statistical  distribution  which 
defines  the  subsurface  heterogeneity,  the  correlation 
distance  of  the  variability,  and  the  mean  and  standard 
deviation  of  the  material  property  being  considered.  These 
parameters  were  chosen  based  on  the  results  of  a  field 
investigation  effort  conducted  in  dry  alluvium. 

A  technique  was  developed  to  construct  spatially  random 
geologic  variability  factors  which  were  Incorporated  into 
the  AFTON  1>D  finite  difference  code  to  simulate  the 
subsurface  inhomogeneities.  This  technique  was  named  the 
Random  Geology  Generator  in  this  study. 


4.2  Implications 

The  results  of  the  calculations  which  incorporated  the 
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spatial  subsurface  variability  Into  the  geologic  material 
model  have  several  Important  Implications  for  the  explosive 
effects  community.  The  Implications  are  as  follows: 


1.  The  calculations  have  shown  that  relatively  small 
percentages  of  variability  (5  percent  standard 
deviation)  result  In  changes  In  the  ground  motion 
waveforms  relative  to  the  homogeneous  calculations, 
binder standing  the  effects  that  spatial  subsurface 
variability  has  on  ground  motion  will  reduce  the 
uncertainty  Involved  when  performing  ground  shook 
oalculations/predlotlons. 

\ 

2.  The  introduction  of  spatial  random  variability  into 
the  elastic  material  caused  broadening  of  peaks,  and  the 
attenuation  rate  to  Increase  relative  to  the  homogeneous 
calculations.  Dynamic  material  model  parameters  are 
often  obtained  by  adjusting  finite  difference  code  Input 
parameters  until  the  calculated  waveforms  match 
experimental  waveforms  In  terms  of  attenuation  rate,  and 
pulse  width.  By  adjusting  certain  parameters,  the 
Increased  attenuation  rate  and  pulse  broadening  might  be 
In  correctly  modeled  as  non-linear  material  effects, 
when  they  can  be  caused  by  spatial  variability  In  a 
linear-elastic  material. 
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3.  The  presence  of  spatial  variability  neans  a  single 
gage  record  at  a  given  range  for  an  experiment  does  not 
provide  a  definitive  view  of  the  ground  motion  field. 
Multiple  gage  records  at  the  same  range  on  different 
azimuths  are  required  to  adequately  define  the  ground 
notion  uncertainty  and  variability. 

4.  Incorporation  of  spatial  variability  into 
calculatlonal  predictions  of  explosive  effects  on 
structures  would  better  define  the  uncertainty  in  the 
predictions.  For  structures  of  large  areal  extent, 
spatial  variability  in  the  ground  motion  field  means 
that  different  portions  of  the  structure  are  subjected 
to  different  loadings. 

5.  Random  spatial  variability  nay  have  a  significant 
influence  on  the  determination  of  explosive  yields  under 
nuclear  testing  treaties  using  near-source  ground  motion 
data.  The  calculations  performed  in  this  study  have 
shown  an  increase  in  the  attenuation  rate  relative  to 
the  homogeneous  cases.  In  addition,  there  is  a 
significant  uncertainty  band  (maximum  band  was  ±  25 
percent)  for  the  velocity  amplitudes.  Again  a  single 
gage  at  a  given  range  may  not  adequately  define  this 
uncertainty.  Near-source  heterogeneity  may  also  have 
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strong  Influences  on  the  distant  seismic  signals 
recorded  from  underground  nuclear  tests  (Gupta,  Lynnes, 
and  Wagner,  1991) . 


±i2 _ Rgggmmflndationi 

since  a  site  characterization  effort  cannot  be  done  In  such 
detail  to  determine  each  and  every  heterogeneity  In  the 
subsurface  material  for  a  test  site,  a  new  site  must  be 
characterized  statistically.  Therefore,  a  methodology  for 
statistical  In  situ  characterization  must  be  developed  to 
account  for  the  spatial  subsurface  variability. 

since  the  geologic  material  behavior  In  explosions  Is  In 
general  non-linear,  the  1-dlmenslonal  modeling  techniques 
developed  In  this  study  should  be  applied  to  an  elastic- 
plastic  material  model.  Therefore,  a  methodology  for 
estimating  the  variability  In  the  material  parameters  used 
in  an  elastic-plastic  model  must  be  developed. 

Following  this  step,  the  modeling  techniques  should  be 
extended  to  a  two-dimensional  environment.  A  two- 
dimensional  random  geology  generator  should  be  developed 
and  Incorporated  Into  a  two-dimensional  ground  shock  code. 
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APPENDIX  A 


RANDOM  GEOLOGY  GENERATOR  PROGRAM 


This  Lb  the  listing  for  the  random  geology  generator  program. 

Thia  program  was  written  for  uae  in  the  PC-MATLAB  aoftware.  It 
requires 

no  formal  coding  language. 

function  X  randvar(diat,N.Mn,Sd,A,dg,s) 

% 

%  The  Gaussian  and  Exponential  correlation  functions  are  defined  in 
%  the  frequency  domain.  The  random  variables  are  generated  through 
%  the  random  number  generator  (RAND)  and  then  fourier  transformed 
%  (FFT) . 


%  INPUT  PARAMETES: 

%  diet  -  Gaussian  (1),  and  Exponential  (2) 

%  N  -  number  of  random  variables  to  be  generated 

%  Hn  -  the  mean 

%  Sd  -  standard  deviation 

%  A  -  number  of  sample  Intervals 

%  dg  >  spatial  sample  interval 

%  (NOTE:  a  -  A  *  dg  -  correlation  distance) 

%  a  "  seed 

%*********i»*****#***ilWf**********ilWf*****i»r***************************** 

%  Correlation  Functions 

for  h-0:N/2; 

Tpi-2.*pl.*h; 
if  (diat  —  1) 

f(h+l)-A.*dg.*sqrt(pi).*exp(-(A*2.*Tpi"2./(N*2.))./4): 
elself  (diet  —  2) 

f(h+l)-(2.*A*dg)/(l+(Tpi/N)*2.*A"2.); 

end 

end 

% 

v-1 

for  h-N/2+l:N-l; 

f  (h+l)-f(N/24-l-v) ; 
v-v+1; 
end 

% 

dlspChit  any  key  to  plot  Theoretical  Transform  of  Autocorrelation') 
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pause 

plot(f) 

% 

%  Writing  output  to  file  called  auto. out 

% 

for  1  -  l:N/2: 

fprlntf( 'auto. out' , '%e  %e\n' ,l,f(l)) 

end 

pause 

t 

%  Random  Number  Generator 

% 

rand(* normal') ; 
rand( 'seed' ,s) 
x-rand(l,N) ; 

% 

%  FFT  on  random  numbers 

% 

x«fft(x) ; 
y-x.*conj (x) ; 

% 

dlspChlt  any  key  to  display  Power  Spectrum  of  Random  Numbers') 
pause 

plot(y) 

pause 

% 

%  Realizations  constructed  by  taking  the  square  root  of  the 
%  power  spectrum  and  multiplying  It  to  the  phase  term  exp(l*phl) 

% 

phl-atan2(lmag<x).real(x)) ; 
x-sqrt<f) .*exp(aqrt(-l) .*<phl)); 
y-x.*conj (x) ; 

t 

dlspChlt  any  key  to  display  Filtered  Power  Spectrum  of  Random 
Numbers ' ) 

pause 

plot(y) 

pause 

% 

%  IFFT  required  to  transform  random  numbers  to  the  space  domain 

% 

x-lfft(x) ; 
x-real(x) ; 

% 

%  Data  Standardization:  Transforming  to  match  the  original 
%  standard  deviation 

% 

x-(x-mjan(x)) .*Sd./std(x<l:N))+Mn; 

% 

dlspChlt  any  key  to  plot  Filtered  Random  Numbers') 
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pause 

plot(x) 

pause 

% 

%  Writing  random  number  output  to  file  called  gaus.dut  for  Gaussian 
%  Correlation  and  exp. out  for  Exponential  Correlation 

t 

for  i-l:N; 

if  (diet  —  1) 

fprlntf ( 'gaus .out* , '%e  %e\n' , l,x(i)) 

fprlntf( 'gausb.out' , 'ta\n' ,x(l)) 
elseif  (dlst  —  2) 

fprlntf ( 'exp. out' , '%e  %e\n' , l,x(i) ) 

fprlntf ( 'expb. out' , '%e\n' ,x(l)) 
and 

end 

pause 

% 

%  Construct  autocorrelation  by  FFT  of  padded  sequence 

% 

xffl-mean<x) ; 
x-(x-xm) ; 
x-fft(x,N*2): 
x-x,*conJ (x) ; 
x-lfft<x) ; 
x-real(x) ; 
x-x./x<l) ; 

% 

for  h-0:N*2-l 
r(h+l)-h.*dg! 

end 

% 

for  ec-0:N*2-l 
rl(ec+l)*ec .*dg; 
if  (diet  —  1) 

ex(ec+l)-exp(-<ec''2.*dg*2.)./((A*dg)'‘2. )) : 
elseif  (dist  —  2) 

ex<ec+l)-exp(- (ec*dg)/(A*dg)) ; 

end 

end 

% 

dispChit  any  key  to  display  Final  Autocorrelation') 

pause 

plot(rl,ex,r,x) 

pause 

%  Writing  theoretical  and  computed  autocorrelation  to  file  final. out 

% 

for  1-1 :N; 

fprlntf (' flnal.out' , '%e  %e  %e\n’ ,r(i) ,ex(l) ,x(i) ) 

end 
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APPENDIX  B 
CHI-SQUARE  TEST 

B.l  Chi-Sauar«  Teat  Deviations 

Th«  most  direct  way  to  test  a  random  number  generator  is  to 
use  It  to  generate  the  desired  random  numbers,  u/s,  and 
examine  the  U/s  statistically  to  determine  how  closely  they 
appear  to  be  uniformly  distributed  between  0  and  1.  The 
chi-square  test,  serial  test,  and  runs  (or  runs-up)  test 
are  a  few  empirical  teste  (Law  and  Kelton,  1982)  .  A 
computer  program  of  the  chi-square  test  of  uniformity  was 
written  to  test  the  function, 

Zi  •  <7*  Zi.^)  (mod  2*’-  -  1)  (8.1) 

The  following  steps  are  reqpilred  to  apply  the  chi-square 
test: 

1.  The  Interval  (0,  1)  Is  divided  Into  k  subintervals 
of  equal  length,  and  Ut,  U2,..>,Ub  are  generated.  (As  a 
rule,  k  should  be  at  least  100,  and  n/k  should  be  at 
least  5.) 

2.  For  j  «  1,  2,...,k,  let  fj  be  the  number  of  the  u/s 
that  are  In  the  jth  subinterval,  and  let 
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(B.2) 


I  S  - 


J  ■  1 


3.  For  large  n,  will  approximately  have  a  chi- 
square  distribution  with  (k  -  1)  df  (degrees  of 
freedom)  under  the  null  hypothesis  that  the  U{'s  are 
IID  U(0,1)  random  variables.  The  null  hypothesis  Is 
rejected  at  level  a  if  x*  >  x\.i,i^»  where  x*n,i<i  is  the 
upper  1-a  critical  point  of  the  chl-scpiare 
distribution  with  k-l  df. 

4.  For  large  values  of  k,  the  approximation 

xi.,.,..  -  (Jt  .  1)  (l  -  .  r,..  (jTiVirn  <»•  »> 

can  be  used,  where  Is  the  upper  1-a  critical 
point  of  the  N(0,1)  distribution. 

The  Input  parameters  used  In  the  chi-square  test  are: 

1.  The  Initial  seed:  Zq  -  12,345,678 

2.  k  -  2*’  -  4096 

3.  n  -  2“  -  65,536 
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The  approximation  for  the  critical  point  at  avO.lO  is: 


*  ik  -  1)  |l 


2  ^  /  2 
9  (ic  -  1)  \9(;f  -  1)  /  j 


X4095,0.9  *  (4095)  |l 


2 

9  (4095) 


4  1.282 


/  -2 
\9  (4095) 


Xli.i-  •  4211.11 


-  4000.8 


Therefore, 

X*  <  xi-i.i-« 

80  we  fail  to  reject  the  NULL  hypothesis  at  the  0.10  level. 
The  random  numbers  generated  do  not  behave  significantly 
different  from  truly  IID  U(o,l)  random  variables,  based  on 
the  chi-sguared  test. 

The  program  written  to  test  the  random  number  generator 
with  the  chi-squared  test  is  Included  in  the  following 
section. 
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CHI-SQUARE  TEST 


This  listing  Is  for  the  chl-squsre  tsst  program.  The  coding  language 
is  Fortran. 

PROGRAM  RANDOM  NUMBER  GENERATOR  TEST 
c  The  chl-aquare  test  la  designed  to  check  whether  the  random 
o  numbers  (U'a)  appear  to  be  uniformly  distributed  between  0  and 
c  1.  [0,  1]  la  divided  Into  sublntervala  (s)  of  equal  length 

c  and  Ul,  U2,  U3,,,.,Un  are  generated, 
c 
c 

c  List  of  symbols: 
c 

c  bm  subinterval  matrix 

c  RNUM  -  generated  random  numbers  (MATLAB) 
c  B  m  number  of  subintervals 
c  ds  size  of  sub  Interval 

c  n  -  number  of  random  numbers 
c  f  -  summation  of  U's 

c  chl2  «  equation  to  calculate  chi- squared 

c 

c 

double  precision  bm(4096, 2) ,RNUH(66000) , f (0:4096) ,d8,chl2 ,n,s 
c 

open(unlt-10,flle-' rout ’ , s  tatus- 'old') 
open  (uni t>l, file- 'output') 
c 

ds  -  1.0/4096.0 
n  -  65536.0 
8  -  4096.0 
o 

c  This  loop  will  generate  the  values  for  matrix  b  -  subintervals 
c 

do  100  1-1,4096 
bm(l,l)  -  ds  'I'  1 
bm(l,2)  -  0.0 
100  continue 
c 

c  These  loops  will  count  the  number  of  U's  In  each  subinterval 
c 

do  10  1  -  1,65536 

read  (10,*)  RNUM(i) 
c 
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if(RNUM<t)  .It.  0.25)  then 
do  20  k-l,I024 

if(RNUM(l)  .It.  bm(k,l>)  then 
bn(k,2)  -  bni<k,2)  +  1.0 
go  to  10 
endif 

20  continue 

else 

i£(RinJM(i)  .It.  0.50)  then 
do  30  k-102S,2048 

if(RNUM(i)  .It.  bD<k,l))  then 
bB(k,2)  -  ba(k,2)  +  1.0 
go  to  10 
endif 

30  continue 

else 

if(RNUM(i)  .It.  0.75)  then 
do  40  k-2049.3072 

lf(RNUM(l)  .It.  bm(k,l))  then 
ba(k,2)-  bB(k,2)  4-  1.0 
go  to  10 
endif 

40  continue 

else 

do  SO  k-3073.4096 

if  (RNUMd)  .It.  bB(k.l>)  then 
ba(k,2)  -  ba(k,2)  4-1.0 
go  to  10 
endif 

SO  continue 

endif 
endif 
endif 

10  continue 
c 

f<0)  -  O.dO 
c 

c  Generation  of  chi- squared  value 
c 

do  90  J-1,4096 

f(j)  -  (bm(J.2)  -  n/a)**2.d0 

f(J)  -  f(j-l)  f(J) 

90  continue 

c 

chi2  -  s/n  ★  f(4096) 
c 

c  Write  fonaat  stateuents  to  generate  output 
c 

write  (1,60) 

60  format  (lx, 'CHI-SQUARE  TEST  RESULTS;',/) 
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write  (1,65)  chi2 

65  format  (lx, 'chl2  -  ',el7.10,//> 
write  (1,66) 

66  format  ( 5x , ' SUBINTERVAL ' , 6x , ' NUMBER  OF  Us ' ) 
write  (1,70)  ( (bm(i. j ) , J-1 , 2) , i-1 ,4096) 

70  format  (2(lx,el7.10)) 

stop 
and 
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APPENDIX  C 


SHARPE'S  CLOSED  FORM  SOLUTION 


Th«  dlaplacanant  In  tha  nadium  corrasponding  to  the 
application  of  unit  function  praesura  to  tha  intar ior 
surfaoa  of  the  cavity  la  aa  followa  (Sharps,  1942): 

u  ■  ^  l-jj*  +  tan"^  v/2) 

+  \/i  I-jI  ainwtj 

where,  u  ■  diaplacenent 

For  diaplacamenta  at  diataneea  more  than  a  few  tinea  the 
radlua,  Equation  G.l  ia  alnplifiad  to  the  following 
axpreaalon 

U  -  Bin  UT  C.2 

2v'2r<3r 

where, 

a  -  cavity  radius 
Po  •  applied  preaaura 
6  ■  ahear  modulus 

r  ■  range  of  Interest 
(i)  ■  2\/2Cp/3a 

Cp  ■  P-wave  velocity 
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T-  -  t-((r-a)/Cp) 

t  «  time  required  to  reach  specified  range 

The  velocity  equation  wae  obtained  by  taking  the  derivative 
of  the  displacement  equation  C.2  (du/dr) 

n 

V  ■  ^  ttcosuT  +  slnut]  C.3 

dt  2^az 

where,  v  *  velocity 

The  shear  modulus  is  obtained  from  the  following  equations: 


C,.  (!)■'■  0.S 

where , 

N  "  oonstalnt  modulus 
V  •  polsson's  ratio 
p  •  soil  density. 

If  Equation  C.5  is  written  as  M  >  Cp^p  and  substituted  into 
Equation  C.4,  the  following  expression  is  obtained 
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cgpd  -  2V) 
2(1  -V) 


C.6 


Equation  c.6  was  alnpllfiad  by  lotting  tho  polsaon'a  ratio, 
V,  aqual  0.25: 

a  •  0.33  pcj  C.7 

A  MATLAB  prograa  waa  writtan  with  Equations  C.2,  C.3,  and 
C.7  to  ganorata  tha  valocity  and  dlaplacanant  wavaforma 
using  Sharps' s  closad  form  solution. 

Tha  valocity  and  dlaplaoamant  wavafoinns  wars  ganaratad  with 
tha  following  Input  paraaatars: 

cavity  radlua  (a)  ■  175  on 

Appllad  proaauro  (Po)  >■  lOa-fio*  dyna/cm^  (10  Bars) 
Danslty  (p)  >  1.643  qfcw? 

P-wavo  valocity  (C^)  >  50800  cxn/sac 
Numbar  of  tlmostapa  -  300 
Tlmastopa  ■  0.0001  aoc 

Ranga  of  Intarast  (r)  ■  200,  500,  900,  and  2000  cm 
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Su2 _ Clggfl<ii  Fgm  Sglutlgn  frogram 


Sharpe's  Closed  Form  Solution  Program 


This  program  was  written  for  use  In  the  PC*MATHlAfi  software.  It 
require 

no  formal  coding  language. 


[unction  X  -  aharpe2a(crad,po,rho,r,tlna,T,c) 


Sharpes  solution  for  step  function  pressure  applied  to  the  wall 
of  a  spherical  cavity. 


Symbols : 

crad 

«■ 

cavity  radius 

po 

- 

applied  pressure 

rho 

- 

density 

r 

range 

tlno 

■i 

time  step  (make  tine  small  enough  to  catch 
waveform) 

T 

#  of  time  steps  to  be  evaluated 

c 

■i 

P  wave  velocity 

G 

shear  modultis  (rigidity)  obtained  from  following 
equation  (G  «  rho*e*c*(l-2v)/(2(l-v)) ,  where  v  Is 
the  polsson's  ratio  of  0.23) 

u 

- 

displacement  output 

udot 

■i 

velocity  output 

NOTE:  any  units  can  be  used  as  long  as  all  are  consistent 


t  Time,  td,  required  for  wave  to  reach  specified  range  r  is  computed 

% 

td-(r-crad) ./c 
Tdl-td/tlnc 
Tdl>round(Tdl) 
om-2*sqrt(2)*c./(3.*crad) ; 

% 

%  Shear  Modulus  Computation 

I 

G  -  O.33*rho*c*2. 

t 

I  Delay  or  travel  time  Is  tacked  onto  front  of  sharpes  waveform 
t 

for  l-l:(Tdl-l); 
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u(i)-0.0: 

udot(i)-0.0; 

end; 

% 

h-O; 

for  i-Tdl:(T+Tdl) ; 
h"h+l : 

t*u«(h-l)*tlnc; 

% 

%  Eveluetlon  of  equation  13  and  Ite  derivative  from  Sharpes  part  1 
i  follows: 

% 

u<i)-(crad.*2.*po./(2*sqrt<2) .*G.*r)).*exp((-tau.*om)/sqrt(2)) .*8in(t 
au*om) ; 

ud-(crad.*2.*po./<2*aqrt<2) .*G.*r)) ; 
der iv-isxp ( ( •  tau^om) .  /sqrt (2 ) ) , *om . '*'oos ( tau'*'om) ; 
doriv-deriv4>((-om./aqrt(2))  .*axp(*tau*om/aqrt(2))  .*8in(tau*ois)) ; 
udo t ( i ) -ud . *da r iv : 
and 
% 

for  h-  Tdl:(Tdl+T) 
tl(h)-(h-l)*tlno; 
and 

% 

dispChit  any  kay  to  plot  displacement*) 

pause 

ploe(ti,u) 

pause 

disp('hit  any  key  to  plot  velocity') 

pause 

plot(ti,udot) 

pause 

% 

%  Convert  units  to  tine  >  usee;  velocity  -  cm/usec;  displacement  cm; 

% 

for  i  -  1:(T+Tdl) ; 

ti(i)  -  tl(i)./1.0e-6; 
udot(i)  -  udot(i) .*1 .Oe-6; 
end; 

% 

diap( 'displacement  plot  •  cm  v..  sec') 

pause 

ploe(tl,u) 

pause 

dispC 'velocity  plot  -  cm/usec  vs  usee') 
pause 

ploc(ti,udoc) 

pause 
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% 

%  writing  output  to  file  called  sharp. out,  col  1  is  time,  2  is 
%  velocity, 

%  3  Is  displacement 

% 

for  1  -  1:(T+Tdl) 

fprintf( 'sharp. out' , '%e  %e  %e\n' ,tl(i) ,udot(i) ,u(l)) 

end 
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NAVAL  CIVIL  ENGINEERING  LABORATORY 
ATTN:  CODE  L64  LORY 

NAVAL  POSTGRADUATE  SCHOOL 
ATTN:  CODE  1424 

NAVAL  RESEARCH  LABORATORY 
ATTN:  0008  2627 

NAVAL  SURFACE  WARFARE  CENTER 
ATTN;  CODE  K42  R  ROBINSON 
AHN;  CODE  K42  8  HUGHES 

NAVAL  WEAPONS  EVALUATION  FACILITY 
ATTN;  CLASSIFIED  LIBRARY 

NEW  LONDON  LABORATORY 
ATTN:  TECH  LIBRARY 

OFFICE  OF  CHIEF  OF  NAVAL  OPERATIONS 
ATTN;  OP03EG 

OFFICE  OF  NAVAL  RESEARCH 
ATTN;  CODE1132SM 

DEPARTMENT  OF  THE  AIR  FORCE 

AFIS/INT 

ATTN:  INT 

AIR  FORCE  TECHNICAL  APPLICATIONS  CTR 
ATTN:  R HANSEN 

AIR  UNIVERSITY  LIBRARY 
ATTN:  AUL-LSE 

PHILLIPS  LABORATORY 

ATTN;  NTES  LTCOL  T  BRETZ 
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UNITED  STATES  STRATEGIC  COMMAND 
ATTN;  J  633 
ATTN:  JIC/ODM 

DEPARTMENT  OF  ENERGY 

DEPARTMENT  OF  ENERGY 

NEVADA  OPERATIONS  OFFICE 
ATTN;  OTIS  D  H  MARTIN 

LAWRENCE  LIVERMORE  NATIONAL  UB 
ATTN:  L-140  R  DONG 
ATTN:  L-200  D  GLENN 
ATTN:  L-203  R  SCHOCK 
ATTN:  L-84G  POMYKAL 
ATTN:  W  CROWLEY 

LOS  ALAMOS  NATIONAL  LABORATORY 
ATTN:  BOB  DEUPREE 
ATTN:  FREDAPP 
ATTN:  R  P  WEAVER 
ATTN:  D  STROTTMAN 
ATTN:  REPORT  LIBRARY 
ATTN:  8  TAYLOR 
ATTN:  R  KIRBY 

MARTIN  MARIETTA  ENERGY  SYSTEMS  INC 
ATTN:  DR  C  V  CHESTER 

SANOIA  NATIONAL  LABORATORIES 
ATTN;  ACHABAID1V9311 
ATTN:  DIV  5214  JS  PHILLIPS 
ATTN:  DIV  9311  LR  HILL 
ATTN:  OR  CARL  W  SMITH  DIV  931 1 
ATTN:  MIKE  FURNISH 
ATTN:  P  YARRINGTON 
ATTN:  TECH  UB  31 41 
ATTN:  9311  CW  SMITH 
ATTN:  9311DQARBIN 

OTHER  GOVERNMENT 

CENTER  FOR  SEISMIC  STUDIES 
ATTN;  L  GRANT 

CENTRAL  INTELLIGENCE  AGENCY 
ATTN;  OSWR/NED 

DEPARTMENT  OF  THE  INTERIOR 
ATTN:  D  RODDY 

FEDERAL  AVIATION  ADMINISTRATION 
ATTN:  L  MOLINAR 

DEPARTMENT  OF  DEFENSE  CONTRACTORS 

AEROSPACE  CORP 

ATTN:  LIBRARY  ACQUISITION 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN:  C  J  HIGGINS 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN;  S  BLOUIN 

APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN;  R  FRANK 


APPLIED  RESEARCH  ASSOCIATES,  INC 
ATTN:  J  L  DRAKE 

APPLIED  THEORY,  INC 
AHN:  J  TRUUO 

BDM  FEDERAL  INC 

ATTN:  E  DORCHAK 
ATTN:  J  STOCKTON 

BOEING  TECHNICAL  &  MANAGEMENT  SVCS,  INC 
ATTN:  ROBERT  M  SCHMIDT 

CALIFORNIA  INSTITUTE  OF  TECHNOLOGY 
ATTN:  T  AHRENS 

CALIFORNIA  RESEARCH  &  TECHNOLOGY,  INC 
ATTN;  J  THOMSEN 
AHN:  K  KREYENHAQEN 

CARPENTER  RESEARCH  CORP 
ATTN;  H  J  CARPENTER 

HORIZONS  TECHNOLOGY,  INC 
ATTN:  B  LEE 

IIT  RESEARCH  INSTITUTE 
ATTN:  A  BUTI 

ATTN:  DOCUMENTS  LIBRARY 
ATTN:  M  JOHNSON 

INSTITUTE  FOR  DEFENSE  ANALYSES 
ATTN;  CLASSIFIED  LIBRARY 

KAMAN  SCIENCES  CORP 
ATTN:  LMENTE 
ATTN:  LIBRARY 

KAMAN  SCIENCES  CORP 
ATTN:  DASIAC 

KAMAN  SCIENCES  CORPORATION 
ATTN:  DASIAC 

KTECHCORP 

ATTN;  E  GAFFNEY 

LACHEL  AND  ASSOCIATES,  INC 
AHN;  J  BECK 

LOCKHEED  MISSILES  &  SPACE  CO,  INC 
ATTN:  PHILIP  UNDERWOOD 

LOCKHEED  MISSILES  &  SPACE  CO,  INC 
ATTN;  TECH  INFO  CTR 

LOGICON  R  &  D  ASSOCIATES 
ATTN;  C  K  B  LEE 
ATTN.  D  SIMONS 
ATTN:  LIBRARY 

LOGICON  R  &  D  ASSOCIATES 
ATTN;  B  KILLIAN 
ATTN;  J  WALTON 
ATTN:  L  GERMAIN 

MAXWELL  LABORATORIES,  INC 
ATTN;  J  MURPHY 
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MCDONNELL  DOUGLAS  CORPORATION 
ATTN:  RHALPRIN 

MISSION  RESEARCH  CORP 

ATTN;  TECH  LIBRARY 

NEW  MEXICO  STATE  UNIVERSITY 
2CYS  ATTN:  A  A  MARTINEZ 
ATTN:  A  K  AYER 
ATTN;  TWARD 
ATTN:  W  MCCARTHY 

PACIFIC-SIERRA  RESEARCH  CORP 
ATTN;  HBRODE 
ATTN:  L  JOHNSON 

S-CUBED 

ATTN;  DR  JL  STEVENS 
ATTN:  DR  K  L  MCLAUGHLIN 
ATTN:  KDPYATTJR 
ATTN:  RUFREN2 
ATTN:  S  PEYTON 

S-CUBED 

ATTN:  L KENNEDY 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  DR  M  MCKAY 
ATTN;  H  PRATT 
ATTN:  H  WILSON 
ATTN:  SMELZER 

ATTN:  TECHNICAL  REPORT  SYSTEM 

SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  WLAYSON 


SCIENCE  APPLICATIONS  INTL  CORP 
ATTN:  TECH  LIBRARY 

SOUTHERN  METHODIST  UNIV 

ATTN:  DR  BRIAN  STUMP 

SRI  INTERNATIONAL 

ATTN;  A  FLORENCE 
ATTN;  DKEOUGH 
ATTN;  DR  B  S  HOLMES 

TECH  REPS.  INC 

ATTN;  FMCMULLAN 

TERRA  TEK,  INC 

ATTN:  C  FELICE 

TITAN  CORPORATION 
AHN;  LIBRARY 
ATTN;  S  SCHUSTER 

TRW  SPACE  &  DEFENSE  SECTOR 
ATTN;  W  WAMPLER 

WEIDLINQER  ASSOC,  INC 
ATTN:  H  LEVINE 

WEIDLINQER  ASSOCIATES,  INC 
AHN:  TDEEVY 

WEIDLINQER  ASSOCIATES,  INC 
ATTN;  J  WRIGHT 
ATTN:  M  BARON 

FOREIGN 

BUNDESAN  STALT  FUR  QEOWISSEN  SCHAFFEN 
ATTN;  S  REAMER 
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